/*
 * Copyright 2009-2019 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
#define BOOST_TEST_MAIN

#define BOOST_TEST_MODULE sigma_test
#include <boost/test/unit_test.hpp>
#include <fstream>
#include <votca/xtp/aobasis.h>
#include <votca/xtp/logger.h>
#include <votca/xtp/orbitals.h>
#include <votca/xtp/ppm.h>
#include <votca/xtp/rpa.h>
#include <votca/xtp/sigma_ppm.h>
#include <votca/xtp/threecenter.h>

using namespace votca::xtp;
using namespace std;

BOOST_AUTO_TEST_SUITE(sigma_test)

BOOST_AUTO_TEST_CASE(sigma_full) {

  ofstream xyzfile("molecule.xyz");
  xyzfile << " 5" << endl;
  xyzfile << " methane" << endl;
  xyzfile << " C            .000000     .000000     .000000" << endl;
  xyzfile << " H            .629118     .629118     .629118" << endl;
  xyzfile << " H           -.629118    -.629118     .629118" << endl;
  xyzfile << " H            .629118    -.629118    -.629118" << endl;
  xyzfile << " H           -.629118     .629118    -.629118" << endl;
  xyzfile.close();

  ofstream basisfile("3-21G.xml");
  basisfile << "<basis name=\"3-21G\">" << endl;
  basisfile << "  <element name=\"H\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"5.447178e+00\">" << endl;
  basisfile << "        <contractions factor=\"1.562850e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"8.245470e-01\">" << endl;
  basisfile << "        <contractions factor=\"9.046910e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.831920e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "  <element name=\"C\">" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"S\">" << endl;
  basisfile << "      <constant decay=\"1.722560e+02\">" << endl;
  basisfile << "        <contractions factor=\"6.176690e-02\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"2.591090e+01\">" << endl;
  basisfile << "        <contractions factor=\"3.587940e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"5.533350e+00\">" << endl;
  basisfile << "        <contractions factor=\"7.007130e-01\" type=\"S\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"3.664980e+00\">" << endl;
  basisfile << "        <contractions factor=\"-3.958970e-01\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"2.364600e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "      <constant decay=\"7.705450e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.215840e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"8.606190e-01\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "    <shell scale=\"1.0\" type=\"SP\">" << endl;
  basisfile << "      <constant decay=\"1.958570e-01\">" << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"S\"/>"
            << endl;
  basisfile << "        <contractions factor=\"1.000000e+00\" type=\"P\"/>"
            << endl;
  basisfile << "      </constant>" << endl;
  basisfile << "    </shell>" << endl;
  basisfile << "  </element>" << endl;
  basisfile << "</basis>" << endl;
  basisfile.close();

  Orbitals orbitals;
  orbitals.QMAtoms().LoadFromFile("molecule.xyz");
  BasisSet basis;
  basis.Load("3-21G.xml");

  AOBasis aobasis;
  aobasis.Fill(basis, orbitals.QMAtoms());

  Eigen::MatrixXd MOs = Eigen::MatrixXd::Zero(17, 17);
  MOs << -0.00761992, -4.69664e-13, 8.35009e-15, -1.15214e-14, -0.0156169,
      -2.23157e-12, 1.52916e-14, 2.10997e-15, 8.21478e-15, 3.18517e-15,
      2.89043e-13, -0.00949189, 1.95787e-12, 1.22168e-14, -2.63092e-15,
      -0.22227, 1.00844, 0.233602, -3.18103e-12, 4.05093e-14, -4.70943e-14,
      0.1578, 4.75897e-11, -1.87447e-13, -1.02418e-14, 6.44484e-14, -2.6602e-14,
      6.5906e-12, -0.281033, -6.67755e-12, 2.70339e-14, -9.78783e-14, -1.94373,
      -0.36629, -1.63678e-13, -0.22745, -0.054851, 0.30351, 3.78688e-11,
      -0.201627, -0.158318, -0.233561, -0.0509347, -0.650424, 0.452606,
      -5.88565e-11, 0.453936, -0.165715, -0.619056, 7.0149e-12, 2.395e-14,
      -4.51653e-14, -0.216509, 0.296975, -0.108582, 3.79159e-11, -0.199301,
      0.283114, -0.0198557, 0.584622, 0.275311, 0.461431, -5.93732e-11,
      0.453057, 0.619523, 0.166374, 7.13235e-12, 2.56811e-14, -9.0903e-14,
      -0.21966, -0.235919, -0.207249, 3.75979e-11, -0.199736, -0.122681,
      0.255585, -0.534902, 0.362837, 0.461224, -5.91028e-11, 0.453245,
      -0.453298, 0.453695, 7.01644e-12, 2.60987e-14, 0.480866, 1.8992e-11,
      -2.56795e-13, 4.14571e-13, 2.2709, 4.78615e-10, -2.39153e-12,
      -2.53852e-13, -2.15605e-13, -2.80359e-13, 7.00137e-12, 0.145171,
      -1.96136e-11, -2.24876e-13, -2.57294e-14, 4.04176, 0.193617, -1.64421e-12,
      -0.182159, -0.0439288, 0.243073, 1.80753e-10, -0.764779, -0.600505,
      -0.885907, 0.0862014, 1.10077, -0.765985, 6.65828e-11, -0.579266,
      0.211468, 0.789976, -1.41532e-11, -1.29659e-13, -1.64105e-12, -0.173397,
      0.23784, -0.0869607, 1.80537e-10, -0.755957, 1.07386, -0.0753135,
      -0.989408, -0.465933, -0.78092, 6.72256e-11, -0.578145, -0.790571,
      -0.212309, -1.42443e-11, -1.31306e-13, -1.63849e-12, -0.17592, -0.188941,
      -0.165981, 1.79403e-10, -0.757606, -0.465334, 0.969444, 0.905262,
      -0.61406, -0.78057, 6.69453e-11, -0.578385, 0.578453, -0.578959,
      -1.40917e-11, -1.31002e-13, 0.129798, -0.274485, 0.00256652, -0.00509635,
      -0.0118465, 0.141392, -0.000497905, -0.000510338, -0.000526798,
      -0.00532572, 0.596595, 0.65313, -0.964582, -0.000361559, -0.000717866,
      -0.195084, 0.0246232, 0.0541331, -0.255228, 0.00238646, -0.0047388,
      -0.88576, 1.68364, -0.00592888, -0.00607692, -9.5047e-05, -0.000960887,
      0.10764, -0.362701, 1.53456, 0.000575205, 0.00114206, -0.793844,
      -0.035336, 0.129798, 0.0863299, -0.0479412, 0.25617, -0.0118465,
      -0.0464689, 0.0750316, 0.110468, -0.0436647, -0.558989, -0.203909,
      0.65313, 0.320785, 0.235387, 0.878697, -0.195084, 0.0246232, 0.0541331,
      0.0802732, -0.0445777, 0.238198, -0.88576, -0.553335, 0.893449, 1.31541,
      -0.00787816, -0.100855, -0.0367902, -0.362701, -0.510338, -0.374479,
      -1.39792, -0.793844, -0.035336, 0.129798, 0.0927742, -0.197727, -0.166347,
      -0.0118465, -0.0473592, 0.0582544, -0.119815, -0.463559, 0.320126,
      -0.196433, 0.65313, 0.321765, 0.643254, -0.642737, -0.195084, 0.0246232,
      0.0541331, 0.0862654, -0.183855, -0.154677, -0.88576, -0.563936, 0.693672,
      -1.42672, -0.0836372, 0.0577585, -0.0354411, -0.362701, -0.511897,
      -1.02335, 1.02253, -0.793844, -0.035336, 0.129798, 0.0953806, 0.243102,
      -0.0847266, -0.0118465, -0.0475639, -0.132788, 0.00985812, 0.507751,
      0.244188, -0.196253, 0.65313, 0.322032, -0.87828, -0.235242, -0.195084,
      0.0246232, 0.0541331, 0.088689, 0.226046, -0.0787824, -0.88576, -0.566373,
      -1.58119, 0.117387, 0.0916104, 0.0440574, -0.0354087, -0.362701,
      -0.512321, 1.39726, 0.374248, -0.793844, -0.035336;

  Eigen::VectorXd mo_energy = Eigen::VectorXd::Zero(17);
  mo_energy << -0.612601, -0.341755, -0.341755, -0.341755, 0.137304, 0.16678,
      0.16678, 0.16678, 0.671592, 0.671592, 0.671592, 0.974255, 1.01205,
      1.01205, 1.01205, 1.64823, 19.4429;
  Logger log;
  TCMatrix_gwbse Mmn{log};
  Mmn.Initialize(aobasis.AOBasisSize(), 0, 16, 0, 16);
  Mmn.Fill(aobasis, aobasis, MOs);

  RPA rpa(log, Mmn);
  rpa.configure(4, 0, 16);
  rpa.setRPAInputEnergies(mo_energy);

  Sigma_PPM sigma = Sigma_PPM(Mmn, rpa);

  Sigma_PPM::options opt;
  opt.homo = 4;
  opt.qpmin = 0;
  opt.qpmax = 16;
  opt.rpamin = 0;
  opt.rpamax = 16;
  sigma.configure(opt);

  Eigen::MatrixXd x = sigma.CalcExchangeMatrix();

  Eigen::MatrixXd x_ref = Eigen::MatrixXd::Zero(17, 17);
  x_ref << -0.898354, 5.71768e-07, -1.37292e-08, 1.98142e-07, -0.149915,
      2.16108e-06, -1.16525e-06, -3.17829e-07, -1.88839e-07, 2.21419e-07,
      -2.98344e-08, -0.0269881, -2.13239e-06, -3.09447e-06, -4.98945e-08,
      0.11249, -0.00337645, 5.71768e-07, -0.690315, 4.39806e-07, -3.6005e-07,
      -5.82737e-08, -0.0266566, -0.000334105, -0.000256846, -0.00109173,
      -0.00375978, 0.132508, 1.1112e-07, -0.0402288, 0.000418149, 0.000685308,
      -1.62313e-08, 4.06964e-09, -1.37292e-08, 4.39806e-07, -0.690316,
      1.49062e-07, 5.35669e-09, 0.000201816, -0.0244124, 0.0107079, -0.13183,
      -0.0138696, -0.00147968, -1.15888e-07, 0.000363237, 0.0401121,
      -0.00316551, 7.90554e-09, 3.08003e-09, 1.98142e-07, -3.6005e-07,
      1.49062e-07, -0.690315, -1.48218e-08, -0.000368948, 0.0107048, 0.0244113,
      -0.0139054, 0.131785, 0.00362461, 3.26075e-08, -0.000716146, -0.00315591,
      -0.0401082, -7.83667e-09, 1.75034e-09, -0.149915, -5.82737e-08,
      5.35669e-09, -1.48218e-08, -0.419561, 2.94693e-08, 1.4477e-08,
      2.02462e-10, 1.3573e-08, -1.31684e-08, 7.55841e-10, 0.111589,
      -1.43489e-08, 1.23778e-08, -9.2956e-09, 0.101203, -0.00966821,
      2.16108e-06, -0.0266566, 0.000201816, -0.000368948, 2.94693e-08,
      -0.162521, -2.50885e-08, 2.37541e-08, -6.18131e-05, -0.000396018,
      0.0286842, 7.29432e-08, 0.0853922, -0.000150722, -0.000326281,
      -1.86634e-07, -1.08635e-09, -1.16525e-06, -0.000334105, -0.0244124,
      0.0107048, 1.4477e-08, -2.50885e-08, -0.162521, -6.437e-08, -0.0249187,
      -0.01421, -0.000249755, -3.68182e-08, -0.000249614, -0.0806557,
      -0.0280463, 8.95564e-08, 1.72272e-09, -3.17829e-07, -0.000256846,
      0.0107079, 0.0244113, 2.02462e-10, 2.37541e-08, -6.437e-08, -0.162521,
      0.0142121, -0.0249168, -0.000313469, -4.22882e-08, -0.00025868, 0.0280471,
      -0.0806557, 1.90873e-08, 8.04122e-10, -1.88839e-07, -0.00109173, -0.13183,
      -0.0139054, 1.3573e-08, -6.18131e-05, -0.0249187, 0.0142121, -0.112664,
      5.15694e-09, 2.67503e-08, -6.82025e-08, -3.58992e-06, 0.00311123,
      -0.000577989, 1.57751e-08, 2.46073e-09, 2.21419e-07, -0.00375978,
      -0.0138696, 0.131785, -1.31684e-08, -0.000396018, -0.01421, -0.0249168,
      5.15694e-09, -0.112664, -3.26454e-09, 2.92286e-08, -3.07662e-05,
      0.000577593, 0.00311108, -4.33364e-09, -2.83729e-09, -2.98344e-08,
      0.132508, -0.00147968, 0.00362461, 7.55841e-10, 0.0286842, -0.000249755,
      -0.000313469, 2.67503e-08, -3.26454e-09, -0.112664, -2.55154e-08,
      0.00316414, 9.12634e-06, 2.94706e-05, 6.29731e-09, -3.74354e-09,
      -0.0269881, 1.1112e-07, -1.15888e-07, 3.26075e-08, 0.111589, 7.29432e-08,
      -3.68182e-08, -4.22882e-08, -6.82025e-08, 2.92286e-08, -2.55154e-08,
      -0.17259, -5.78498e-08, -9.45104e-08, 2.34917e-08, -0.0427377, 0.00896338,
      -2.13239e-06, -0.0402288, 0.000363237, -0.000716146, -1.43489e-08,
      0.0853922, -0.000249614, -0.00025868, -3.58992e-06, -3.07662e-05,
      0.00316414, -5.78498e-08, -0.131712, -1.06678e-08, 2.79713e-08,
      1.92221e-07, -7.78508e-10, -3.09447e-06, 0.000418149, 0.0401121,
      -0.00315591, 1.23778e-08, -0.000150722, -0.0806557, 0.0280471, 0.00311123,
      0.000577593, 9.12634e-06, -9.45104e-08, -1.06678e-08, -0.131712,
      2.65503e-07, 2.62621e-07, 2.69838e-10, -4.98945e-08, 0.000685308,
      -0.00316551, -0.0401082, -9.2956e-09, -0.000326281, -0.0280463,
      -0.0806557, -0.000577989, 0.00311108, 2.94706e-05, 2.34917e-08,
      2.79713e-08, 2.65503e-07, -0.131712, 5.71293e-09, -9.80113e-10, 0.11249,
      -1.62313e-08, 7.90554e-09, -7.83667e-09, 0.101203, -1.86634e-07,
      8.95564e-08, 1.90873e-08, 1.57751e-08, -4.33364e-09, 6.29731e-09,
      -0.0427377, 1.92221e-07, 2.62621e-07, 5.71293e-09, -0.10267, 0.0171175,
      -0.00337645, 4.06964e-09, 3.08003e-09, 1.75034e-09, -0.00966821,
      -1.08635e-09, 1.72272e-09, 8.04122e-10, 2.46073e-09, -2.83729e-09,
      -3.74354e-09, 0.00896338, -7.78508e-10, 2.69838e-10, -9.80113e-10,
      0.0171175, -0.0285864;

  bool check_x = x_ref.isApprox(x, 1e-5);
  if (!check_x) {
    cout << "Sigma X" << endl;
    cout << x << endl;
    cout << "Sigma X ref" << endl;
    cout << x_ref << endl;
  }
  BOOST_CHECK_EQUAL(check_x, true);

  sigma.PrepareScreening();
  Eigen::MatrixXd c = sigma.CalcCorrelationOffDiag(mo_energy);
  c.diagonal() = sigma.CalcCorrelationDiag(mo_energy);

  Eigen::MatrixXd c_ref = Eigen::MatrixXd::Zero(17, 17);
  c_ref << 2.0054, -1.80868e-05, 1.13795e-05, -5.53939e-06, -0.0993838,
      1.16048e-05, -8.45082e-06, 3.3849e-07, 4.77622e-06, 2.21006e-06,
      7.61997e-06, 0.265943, -7.04564e-06, -2.62336e-06, -2.09717e-06,
      -0.114711, -0.0189675, -1.80868e-05, 0.0626847, -1.2439e-07, 3.54057e-08,
      -2.44692e-08, 0.000233711, 2.89302e-06, 2.2628e-06, -6.15508e-05,
      -0.000211978, 0.00747013, 1.27237e-06, 0.0364669, -0.00037963,
      -0.00062098, -4.31069e-07, -1.3664e-09, 1.13795e-05, -1.2439e-07,
      0.0626851, 9.39537e-08, -7.87497e-09, -1.82098e-06, 0.000214185,
      -9.39874e-05, -0.00743185, -0.000781895, -8.3411e-05, -6.09584e-07,
      -0.000329304, -0.0363593, 0.00286851, 3.60967e-07, 2.76784e-09,
      -5.53939e-06, 3.54057e-08, 9.39537e-08, 0.0626852, -5.58515e-09,
      3.26137e-06, -9.39016e-05, -0.000214236, -0.000783919, 0.00742934,
      0.000204376, 5.15012e-07, 0.000648951, 0.00286262, 0.0363553,
      -4.71173e-08, 2.7192e-10, -0.0993838, -2.44692e-08, -7.87497e-09,
      -5.58515e-09, -0.0451817, -2.18629e-07, 1.28638e-07, 2.19834e-08,
      -5.8295e-08, -5.75912e-09, -5.27221e-08, 0.11602, 4.80024e-07,
      2.90399e-07, 1.22297e-07, 0.021668, -0.00771886, 1.16048e-05, 0.000233711,
      -1.82098e-06, 3.26137e-06, -2.18629e-07, -0.010061, 9.87618e-08,
      1.40981e-08, -1.50093e-05, -9.5958e-05, 0.00695142, 1.4892e-06, 0.0633079,
      -0.000111616, -0.000241882, -1.20347e-07, 3.7324e-08, -8.45082e-06,
      2.89302e-06, 0.000214185, -9.39016e-05, 1.28638e-07, 9.87618e-08,
      -0.0100613, -6.99123e-08, -0.00603885, -0.00344371, -6.05334e-05,
      -8.1245e-07, -0.000185172, -0.0597966, -0.020793, 1.05236e-07,
      -2.4424e-08, 3.3849e-07, 2.2628e-06, -9.39874e-05, -0.000214236,
      2.19834e-08, 1.40981e-08, -6.99123e-08, -0.0100616, 0.00344415,
      -0.00603842, -7.59713e-05, 4.63477e-08, -0.000191725, 0.0207936,
      -0.059797, -9.88694e-08, 2.22078e-09, 4.77622e-06, -6.15508e-05,
      -0.00743185, -0.000783919, -5.8295e-08, -1.50093e-05, -0.00603885,
      0.00344415, -0.0452211, 3.04907e-08, -2.76997e-09, -6.12261e-07,
      2.50221e-05, -0.0220489, 0.00409566, 5.75757e-07, -1.10819e-08,
      2.21006e-06, -0.000211978, -0.000781895, 0.00742934, -5.75912e-09,
      -9.5958e-05, -0.00344371, -0.00603842, 3.04907e-08, -0.0452213,
      -4.12523e-08, -4.56508e-07, 0.000218076, -0.0040953, -0.0220481,
      2.8037e-07, -4.13616e-09, 7.61997e-06, 0.00747013, -8.3411e-05,
      0.000204376, -5.27221e-08, 0.00695142, -6.05334e-05, -7.59713e-05,
      -2.76997e-09, -4.12523e-08, -0.0452211, -1.11999e-06, -0.0224251,
      -6.44066e-05, -0.000209975, 6.98266e-07, -1.46997e-08, 0.265943,
      1.27237e-06, -6.09584e-07, 5.15012e-07, 0.11602, 1.4892e-06, -8.1245e-07,
      4.63477e-08, -6.12261e-07, -4.56508e-07, -1.11999e-06, 0.154173,
      2.17234e-06, 8.89609e-07, 7.68307e-07, -0.160993, 0.0217415, -7.04564e-06,
      0.0364669, -0.000329304, 0.000648951, 4.80024e-07, 0.0633079,
      -0.000185172, -0.000191725, 2.50221e-05, 0.000218076, -0.0224251,
      2.17234e-06, -0.11857, -2.18742e-07, 8.6041e-08, -2.29157e-07,
      -8.67786e-08, -2.62336e-06, -0.00037963, -0.0363593, 0.00286262,
      2.90399e-07, -0.000111616, -0.0597966, 0.0207936, -0.0220489, -0.0040953,
      -6.44066e-05, 8.89609e-07, -2.18742e-07, -0.11857, 1.81131e-07,
      -5.98803e-08, -9.79642e-08, -2.09717e-06, -0.00062098, 0.00286851,
      0.0363553, 1.22297e-07, -0.000241882, -0.020793, -0.059797, 0.00409566,
      -0.0220481, -0.000209975, 7.68307e-07, 8.6041e-08, 1.81131e-07, -0.118571,
      -2.25477e-07, 3.10441e-09, -0.114711, -4.31069e-07, 3.60967e-07,
      -4.71173e-08, 0.021668, -1.20347e-07, 1.05236e-07, -9.88694e-08,
      5.75757e-07, 2.8037e-07, 6.98266e-07, -0.160993, -2.29157e-07,
      -5.98803e-08, -2.25477e-07, -0.0607528, 0.0159565, -0.0189675,
      -1.3664e-09, 2.76784e-09, 2.7192e-10, -0.00771886, 3.7324e-08,
      -2.4424e-08, 2.22078e-09, -1.10819e-08, -4.13616e-09, -1.46997e-08,
      0.0217415, -8.67786e-08, -9.79642e-08, 3.10441e-09, 0.0159565, -0.408479;

  bool check_c_diag = c.diagonal().isApprox(c_ref.diagonal(), 1e-5);
  if (!check_c_diag) {
    cout << "Sigma C" << endl;
    cout << c.diagonal() << endl;
    cout << "Sigma C ref" << endl;
    cout << c_ref.diagonal() << endl;
  }
  BOOST_CHECK_EQUAL(check_c_diag, true);

  bool check_c = c.isApprox(c_ref, 1e-5);
  if (!check_c) {
    cout << "Sigma C" << endl;
    cout << c << endl;
    cout << "Sigma C ref" << endl;
    cout << c_ref << endl;
  }
  BOOST_CHECK_EQUAL(check_c, true);
}

BOOST_AUTO_TEST_SUITE_END()
